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Executive  Summary 

The  equilibrium  wall  model  is  implemented  in  a  nodal  finite  element  flow  solver  JENRE  developed  at  the  Naval 
Research  Laboratory.  The  Crocco-Busemann  relation  is  used  to  account  for  the  compressibility.  In  this  wall-model 
implementation,  the  first  cell  adjacent  to  the  wall  is  used  to  estimate  the  shear  stress  on  the  wall.  The  no-slip 
adiabatic  boundary  condition  is  applied  to  the  inviscid  and  viscous  fluxes  on  the  wall  to  satisfy  the  surface  physical 
condition,  but  a  non-zero  surface  tangential  velocity  is  used  in  the  calculation  of  the  volume  and  surface  integrals  in 
a  cell  adjacent  to  the  wall.  This  is  because  using  a  zero-surface  velocity  will  grossly  underestimate  these  integrals 
due  to  the  linear  basis  function  used  in  JENRE.  This  implementation  is  validated  in  a  subsonic  boundary-layer  flow 
over  a  flat  plate  and  a  supersonic  flow  in  a  converging  and  diverging  nozzle  used  frequently  in  our  jet  noise 
simulations.  Skin  frictions,  velocity  profiles  and  turbulence  quantities  predicted  by  the  current  wall-model 
implementation  agree  well  with  available  experimental  data  and  theoretical  models.  The  grid  convergence  is 
excellent.  Grid  sizes  much  larger  than  those  recommended  in  other  wall-model  implementations  can  be  used.  The 
current  wall-model  implementation  has  not  encountered  the  numerical  problem  associated  with  using  the  first  cell  as 
reported  in  other  wall-model  implementations.  The  volume  and  surface  integrals  based  on  the  non-zero  surface 
velocity  in  a  cell  adjacent  to  the  wall  show  a  good  agreement  with  those  derived  from  the  equilibrium  boundary- 
layer  velocity  profile  and  the  density  profile  based  on  the  Crocco-Busemann  relation. 
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Equilibrium  Wall-Model  Implementation  in  a  Nodal  Finite  Element 
Flow  Solver  JENRE  for  Large-Eddy  Simulations 

I.  Introduction 

Using  large-eddy  simulations(LES)  to  resolve  wall-bounded  flows  at  high  Reynolds  numbers  is  computationally 
prohibitive  due  to  the  limitations  of  the  numerical  methods  and  computational  resources.  The  major  bottleneck  is  to 
resolve  the  inner  layer  in  the  boundary  layer  near  the  wall.  It  is  a  region  where  most  of  the  turbulence  production 
takes  place  and  turbulent  eddies  are  created.  According  to  Piomelli  and  Balaras  [1],  99%  of  the  grid  points  are  used 
to  resolve  an  inner  layer  which  thickness  is  only  10%  of  the  boundary  layer  at  a  Reynold  number  of  106.  The 
number  of  the  grid  points  required  in  the  outer  layer  increases  with  the  Reynolds  number  as  Reo  s ,  but  the  number  of 
the  grid  points  in  the  inner  layer  increases  as  Re2A  [2].  In  addition,  this  very  fine  grid  resolution  imposes  a  severe 
constraint  on  the  time-step  size.  To  make  LES  applicable  to  wall-bounded  flows  at  high  Reynolds  numbers,  attempt 
to  bypass  the  inner  layer  and  model  its  effects  in  a  global  sense  have  been  used  since  the  pioneering  work  by 
Deardorff  [1]  and  Schumann  [4].  Three  classes  of  modeling  the  inner  layer  have  been  proposed  and  developed  [1]- 
[5].  The  least  expensive  method  is  the  equilibrium  wall-model  approach,  which  assumes  the  existence  of  an 
equilibrium  layer  in  which  the  shear  stress  is  constant.  This  assumption  results  in  the  existence  of  a  logarithmic  layer 
that  can  be  used  to  compute  the  wall  stress  from  the  velocity  in  the  outer  layer.  Other  two  classes  are  the  zonal 
method  and  the  hybrid  method.  Both  of  them  require  the  computation  of  separate,  but  simplified,  flow  equations  in 
the  inner  layer.  More  information  regarding  the  wall-modeling  effort  can  be  found  in  reviews  by  Piomelli  and 
Balaras  [1],  Piomelli  [2]  and  Spalart  [5].  Since  the  wall-model  method  models  the  effect  of  the  inner  layer  on  the  out 
layer,  it  is  possible  to  introduce  some  numerical  errors.  It  has  been  reported  that  the  wall-model  method  often 
encounters  “log-layer  mismatch”  and  10%- 15%  error  in  skin-friction  prediction.  These  types  of  numerical 
discrepancies  can  be  observed  in  all  three  wall-model  implementations.  In  the  equilibrium  wall-model  approach, 
large  numerical  discrepancies  are  often  observed  when  the  first  grid  point  above  the  wall  is  used  to  estimate  the  wall 
shear  stress.  It  has  been  recommended  to  use  grid  points  further  above  the  wall  to  improve  the  prediction  [6]. 

Our  in-house  LES  Solver  JENRE  (Jet  Engine  Noise  Reduction)  uses  a  nodal  finite-element  method  built  upon 
unstructured  meshes.  Since  implementing  the  zonal  and  hybrid  methods  requires  significant  changes  to  the  solver, 
the  equilibrium  wall-model  method  is  used  in  the  current  wall-model  implementation.  It  should  be  mentioned  that 
the  reported  numerical  issues  associated  with  the  wall-model  method  come  from  simulations  using  either  finite- 
difference  or  finite- volume  solvers.  The  implementation  in  finite-element  flow  solvers  for  large-eddy  simulations 
has  rarely  been  reported.  In  this  work,  we  first  examine  the  performance  of  the  equilibrium  wall-mode  method  in  a 
subsonic  boundary-layer  flow  over  a  flat  plate  at  Moo=  0.9.  The  results  are  compared  with  available  experimental 
data.  Since  JENRE  has  been  used  frequently  in  supersonic  jet  noise  simulations  [7]-[  1 0],  this  wall-model  method  is 
also  tested  in  a  supersonic  flow  inside  a  converging  and  diverging  (C-D)  nozzle. 


II.  Large-Eddy  Simulation  Methodology  and  Equilibrium  Wall  Model 

The  compressible  Navier- Stokes  equations  can  be  written  as 

g+V-(F“(<?)-F”  (<?))  =  0,  (1) 

where  Q  is  the  vector  of  the  conservative  variables,  Fa  is  the  inviscid  flux  vector  and  Fvi s  the  viscous  flux  vector. 
The  definitions  of  those  terms  can  be  found  in  any  text  books  of  computational  fluid  dynamics.  Eq.  (1)  is  solved  by 
the  Taylor-Galerkin  nodal  finite  element  method  implemented  in  JENRE.  Taylor-Galerkin  finite  element  method 
maintains  the  second  order  spatial  accuracy  even  on  irregular,  unstructured  tetrahedral  meshes,  which  is  often 
difficult  to  achieve  using  finite  volume  methods.  In  addition,  the  nodal  method  requires  less  computational  memory 
than  the  cell-centered  method  because  it  requires  a  smaller  number  of  degrees  of  freedom  in  tetrahedral  meshes.  In 
the  current  version  of  JENRE,  the  linear  basis  function  is  used  for  the  spatial  discretization,  and  the  finite  element 
flux  corrected  transport  (FEM-FCT)  method  is  used  for  the  flux  integration  and  limiting  calculations  [1 1][12].  The 
multi-dimensional  FEM-FCT  flux  limiter  provides  an  implicit  subgrid  stress  model,  and  the  methodology  is  in  the 
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category  of  the  Monotonically  Integrated  Large-Eddy  Simulation  (MILES)  approach.  The  tetrahedral  mesh  is  used 
because  of  its  simplicity  and  the  accurate  representation  of  complex  geometries. 

The  velocity  profile  in  an  equilibrium  boundary  layer  has  two  distinct  layers  as  shown  in  Figure  1.  The  velocity  near 
the  wall  can  be  modeled  as  a  linear  function  of  the  distance  to  the  wall,  but  slightly  away  from  the  wall,  the  velocity 
distribution  follows  a  logarithmic  function.  The  equilibrium  wall-model  method  uses  the  information  in  the 
logarithmic  layer  to  model  the  effect  of  the  inner  layer  on  the  out  layer.  In  this  approach,  the  viscous  shear  stresses 
on  the  wall  surface  are  not  computed  from  velocity  gradients  near  the  wall.  They  are  estimated  from  the  logarithmic 
velocity  profile  [13]: 


u+  =^-  =  i/n(y+)  +  F,  (2) 

where  k  =  0.41  and  B  =  5.  uT  is  the  friction  velocity  and  ux  =  V tw / pw  .  The  parameters  tw  and  pw  are  the  shear 
stress  and  density  on  the  wall.  In  addition,  u+and  y+  are  the  velocity  and  distance  expressed  in  wall  units,  where 
u  =  u/uT  and  y+  =  pwuTy/pw.  The  parameter  pw  is  the  viscosity  at  the  wall.  The  wall  shear  stress  tw  is  computed 
from  Eq.  (2)  based  on  the  tangential  velocity  U  of  a  given  grid  point  located  in  the  logarithmic  layer.  Because  the 
velocity  and  the  heat  flux  on  the  wall  are  zero  on  an  adiabatic  wall  surface,  the  contribution  of  surface  fluxes  to  the 
flow  is  solely  represented  by  the  wall  shear  stress  in  this  equilibrium  wall-model  method.  Eq.  (2)  is  derived  from  an 
incompressible  velocity  profile,  and  the  density  is  assumed  as  a  constant.  The  compressibility  effect  is  thus  taken 
into  account  by  using  the  van  Driest  effective  velocity  formulation  [13].  The  velocity  U  in  Eq.  (2)  is  replaced  by  the 
van  Driest  effective  velocity  Ueq ,  which  can  be  expressed  as 

Ueq  =  JJC^rZSin-1  -=z=  (3) 

at  an  adiabatic  wall  condition.  Taw  is  the  adiabatic  temperature  on  the  wall,  and  Cp  is  the  specific  heat  capacity.  The 
first  point  above  the  wall  would  be  a  natural  choice  for  the  estimation  of  the  wall  shear  stress.  But  as  mentioned  in 
the  previous  section  that  large  numerical  discrepancies  associated  with  using  the  first  point  have  been  reported  in 
many  wall-model  implementations.  It  has  been  suggested  to  use  points  located  further  above  the  surface,  for 
example,  the  third  or  even  the  fifth  point  from  the  wall  surface  [6]  [15]  [16].  The  implementation  of  this  method  is 
straightforward  in  solvers  using  structured  meshes,  but  it  could  cause  severe  computation  penalties  in  solvers  using 
unstructured  meshes  coupled  with  parallel  computing  methods,  for  example  JENRE.  Thus,  in  the  current  wall-model 
implementation,  the  first  cell  near  the  wall  surface  is  used  to  estimate  the  wall  shear  stress. 

The  no-slip  adiabatic  surface  condition  is  applied  to  the  inviscid  and  viscous  fluxes  on  the  wall  to  satisfy  the  surface 
physic  condition.  Since  the  current  version  of  JENRE  uses  the  linear  basis  function,  using  a  zero-surface  velocity 
would  grossly  underestimate  the  volume  and  surface  integrals  in  a  cell  adjacent  to  the  wall,  for  example,  the 
integrals  of  mass,  momentum  and  kinetic  energy.  To  avoid  this  underestimation,  a  non-zero  tangential  surface 
velocity  is  used  in  those  integrals.  This  tangential  surface  velocity  is  updated  along  with  the  flow  field  at  interior 
nodes.  It  will  be  shown  in  the  following  section  that  the  volume  and  surface  integrals  using  a  non-zero  tangential 
surface  velocity  have  a  good  agreement  with  integrals  derived  from  the  equilibrium  boundary-layer  velocity  profile 
(shown  in  Figure  1)  and  the  density  profile  based  on  the  Crocco-Busemann  relation.  Of  course,  it  would  be  a  more 
ideal  option  to  use  the  equilibrium  boundary-layer  velocity  profile  and  the  Crocco-Busemann  relation  for  the 
velocity  and  density  basis  functions  near  the  wall.  However,  this  would  require  significant  changes  to  the  existing 
code.  It  could  be  considered  in  the  future  development. 


III.  A  Boundary-Layer  Flow  At  M ^  0.9  on  a  Flat  Plate 

Figure  2  shows  the  mesh  distribution  for  this  boundary- layer  flow  over  a  flat  plate.  Inflow  turbulences  are  generated 
by  surface  roughness  presented  as  surface  bumps  in  the  entrance  region.  The  mesh  inside  the  boundary  layer  is 
similar  to  a  structured  mesh  except  that  each  hexahedral  cell  is  divided  into  several  tetrahedral  cells.  This  type  of 
boundary-layer  mesh  makes  it  easier  to  generate  cells  with  large  aspect  ratios,  which  may  reduce  the  number  of  cells 
used  in  boundary  layers  when  length  scales  of  turbulence  structures  are  significantly  different  in  each  direction. 
Another  benefit  of  using  this  type  of  mesh  is  in  the  post  processing,  because  the  location  of  each  node  is  predefined 
and  no  interpolation  is  needed  in  the  data  processing. 
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A.  Using  cubic  cells 

Cubic  cells  where  grid  sizes  are  the  same  in  all  three  directions  are  first  tested.  Table  1  shows  mesh  parameters  of 
four  grid  resolutions  at  the  axial  location  of  Ree  =  Poo^oo^/Moo  =  31000,  where  6  is  the  momentum  thickness. 
Because  the  finest  mesh  requires  a  much  larger  number  of  grid  points,  the  computational  domain  of  the  finest  grid 
resolution  extends  only  as  far  as  Ree  =  22000.  Thus,  the  mesh  parameters  of  this  grid  resolution  shown  in  Table  1 
are  estimated  values.  Figure  3  shows  an  instantaneous  temperature  distribution  using  the  second  coarsest  mesh 
(mesh_II)  shown  in  Table  1.  Disturbances  are  generated  when  the  flow  passes  the  surface  roughness  in  the  entrance 
region  and  the  flow  becomes  turbulent  some  distance  downstream.  The  axial  length  from  the  surface  roughness  area 
to  the  end  of  the  plate  is  1.4  m  in  the  three  coarser  simulations.  The  corresponding  ReQ  is  roughly  44000. 


7.  Comparison  with  the  measurement  data 

Figure  4  shows  streamwise  skin- friction  evolvement  compared  with  an  incompressible  empirical  relation  [18]  and 
incompressible  measurement  data  [19].  To  compare  the  present  compressible  flow  results  with  incompressible 
skin  frictions,  the  van  Driest  II  transformation  [6] [13] [14]  is  used  to  transform  the  compressible  Reynolds  number 
and  skin  friction  into  equivalent  incompressible  values, 


f,inc 


Tw/Tcc 


■Cf,  a  =  V 1  —  Toa/Tw  ,  Ree,inc  =  ^ Ree , 

Mw 


(4) 


where  the  skin  friction  Cf  =  2Tw/(p00f/(^).  The  four  simulations  present  a  very  good  grid  convergence  where  the 
boundary  layer  is  established.  The  magnitude  and  axial  slope  agree  very  well  with  those  shown  in  measurement 
data.  Figure  5  shows  the  comparison  of  streamwise  velocity  distributions.  The  distributions  of  the  three  coarser 
meshes  are  results  at  the  location  of  Ree  inc  =  31000,  where  the  boundary-layer  thickness  is  roughly  0.02  m.  The 
result  of  the  finest  mesh  is  at  the  end  of  its  computational  domain  at  Ree  =  22000.  Again,  there  is  a  very  good  grid 
convergence  in  the  logarithmic  layer,  and  the  predictions  show  a  good  agreement  with  the  measurement  data  and  the 
logarithmic  function.  Figure  6  compares  the  turbulent  intensities  and  the  Reynolds  shear  stress.  The  numerical 
predictions  agree  reasonable  well  with  the  measurement  data.  As  the  grid  resolution  increases,  more  turbulence 
energy  is  captured  and  the  peak  location  shifts  further  towards  the  wall.  The  differences  shown  near  the  boundary- 
layer  edge  are  very  likely  introduced  by  the  differences  in  the  ambient  turbulence  levels.  There  are  ambient 
turbulences  in  measurement  data,  but  the  ambient  turbulences  are  not  included  in  the  numerical  simulations. 


2.  Shape  factor  calculation 

The  shape  factor  is  the  ratio  between  the  displacement  thickness  and  the  momentum  thickness.  It  is  an  important 
factor  to  judge  if  the  boundary-layer  flow  is  laminar  or  turbulent.  Laminar  boundary  layers  have  shape  factors 
around  2.5,  but  turbulent  boundary  layers  have  lower  shape  factors  between  1.3  and  1.8  [20].  Since  the  first  cell  size 
above  the  wall  is  not  small  in  the  wall-model  approach,  how  to  compute  the  mass  and  momentum  integrals  inside 
the  first  cell  can  significantly  affect  the  shape  factor.  Two  approaches  are  used  in  the  shape-factor  calculation.  In  the 
first  approach,  the  velocity  and  density  are  approximated  by  linear  functions  of  the  normal  distance  in  the  first  cell, 
but  in  the  second  approach,  the  velocity  profile  of  the  equilibrium  boundary  layer  shown  in  Figure  1  and  the  Crocco- 
Busemann  relation  are  assumed.  Further  above  the  wall,  the  velocity  and  density  are  approximated  by  linear 
functions  inside  each  cell  to  be  consistent  with  the  methodology  used  in  JENRE.  Figure  7  presents  the  shape-factor 
distributions  computed  by  these  two  approaches.  If  a  linear  approximation  is  used  inside  the  first  cell,  the  shape 
factors  shown  in  Figure  7  (a)  decrease  in  the  axial  direction,  similar  to  the  general  trend  observed  in  measurement 
data.  The  shape  factors  have  large  magnitudes,  and  they  decrease  as  the  grid  resolution  increases.  The  shape  factors 
predicted  by  the  two  finer  meshes  are  below  1.8  at  some  downstream  locations.  If  judged  by  the  results  shown  in 
Figure  7  (a),  the  boundary- layer  flows  predicted  by  the  two  coarser  meshes  are  not  yet  turbulent.  But  if  the  velocity 
profile  of  the  equilibrium  boundary  layer  and  the  density  profile  based  on  the  Crocco-Busemann  relation  are  used  in 
the  first  cell,  the  shape  factor  is  much  less  sensitive  to  the  grid  resolution,  as  shown  in  Figure  7  (b).  The  axial 
variation  is  also  small,  and  the  shape  factor  is  around  1.6  at  all  four  grid  resolutions.  It  would  be  interesting  to  see 
what  the  shape  factor  is  if  the  mass  and  momentum  integrations  are  computed  based  on  the  velocity  profile  of  the 
equilibrium  boundary  layer  and  the  Crocco-Busemann  relation  in  the  entire  boundary  layer.  These  integrations 
require  the  knowledge  of  the  outer  edge  of  the  logarithmic  layer.  The  distance  at  the  outer  edge  of  the  logarithmic 
layer  is  between  5000  and  15000  wall  units  in  the  axial  range  shown  in  Figure  7,  and  the  corresponding  shape  factor 
is  between  1.59  and  1.56,  respectively.  These  values  are  very  close  to  the  shape  factors  shown  in  Figure  7  (b).  It 
should  be  mentioned  that  these  integrations  do  not  include  the  mass  and  momentum  deficits  beyond  the  logarithmic 
layer.  However,  because  the  contribution  of  the  region  beyond  the  logarithmic  layer  should  be  small  in  an 
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equilibrium  turbulent  boundary  layer,  it  is  not  expected  that  the  shape  factor  of  an  equilibrium  turbulent  boundary 
layer  deviates  significantly  from  the  values  shown  above.  It  can  be  thus  concluded  that  the  boundary  layers  in  these 
four  simulations  are  turbulent  in  the  axial  range  shown  in  Figure  7.  This  also  brings  up  an  issue  associated  with  the 
shape-factor  calculation.  It  is  believed  that  the  linear  approximation  between  two  neighboring  vertical  points  has 
been  used  to  compute  the  mass  and  momentum  integrations  in  measurements.  This  may  render  caution  in  using 
measured  shape-factor  data.  The  physical  locations  of  the  measured  points,  especially  the  point  adjacent  to  the  wall, 
should  be  examined  to  see  if  they  are  appropriate  for  a  linear  approximation  between  two  neighboring  points. 

3.  Benefits  of  using  a  non-zero  surface  tangential  velocity 

As  mentioned  above,  a  non-zero  tangential  surface  velocity  is  used  in  the  volume  and  surface  integrations  in  cells 
adjacent  to  the  wall.  Those  integrations  can  be  expressed  as  the  mass  integral  f  p  dy,  the  momentum  integral 
f  pu  dy,  and  the  kinetic  energy  integral  J  pu2  dy  in  the  first  cell,  where  y  is  the  direction  normal  to  the  wall.  It  is 
worthwhile  to  compare  those  integrals  obtained  from  simulations  based  on  a  non-zero  surface  velocity  and  those 
derived  from  the  velocity  profile  of  an  equilibrium  turbulent  boundary  layer  and  the  density  profile  derived  from  the 
Crocco-Busemann  relation.  Figure  8  and  Figure  9  show  comparisons  of  f  pu  dy  and  f  pu2  dy  in  the  first  cell. 
Predictions  using  the  four  grid  resolutions  are  presented,  and  the  results  based  on  the  linear  shape  function  and  the 
zero-surface  velocity  are  also  included.  There  is  an  excellent  grid  convergence  and  the  axial  variations  are  very 
small.  The  momentum  integrals  are  only  around  6%  larger  and  the  kinetic-energy  integrals  roughly  10%  larger  than 
those  derived  from  the  velocity  profile  of  the  equilibrium  turbulent  boundary  layer  and  the  density  profile  base  on 
the  Crocco-Busemann  relation.  This  comparison  should  be  considered  very  good  since  boundary- layer  flows  are  not 
expected  to  follow  those  analytical  functions  exactly.  The  difference  shown  in  the  mass  integrals  f  p  dy  is  even 
smaller  (not  shown),  which  is  less  than  1%.  On  the  other  hand,  if  the  zero-surface  velocity  is  used,  the  momentum 
and  kinetic-energy  integrals  based  on  the  linear  basis  function  are  grossly  underestimated.  These  comparisons 
clearly  demonstrate  the  effectiveness  of  using  a  non-zero  surface  velocity  in  the  wall-model  implementation  in  flow 
solvers  using  the  linear  basis  function. 

B.  Using  cells  of  large  aspect  ratios 

It  is  well  known  that  turbulent  structures  near  the  wall  are  elongated  in  the  streamwise  direction  in  subsonic 
boundary  layers  over  a  flat-plate.  Thus,  the  mesh  type  with  a  larger  streamwise  grid  size  is  expected  to  be 
appropriate.  This  mesh  type  is  tested  and  the  streamwise  grid  size  is  fixed  at  a  large  value  of  Ax+  =  370.  The  grid 
sizes  (Ay  =  Az )  in  other  two  directions  vary.  Results  of  three  aspect  ratios  (Ax+/  Ay+),  1.0,  1.5  and  2.25,  are  shown 
in  Figure  10.  It  can  be  seen  that  there  is  a  good  grid  convergence  and  predictions  agree  well  with  both  the 
measurement  data  and  the  logarithmic  law.  But  if  a  very  large  spanwise  mesh  size  is  also  used  that  A z+  =  Ax+  =  370, 
Figure  1 1  shows  that  the  grid  convergence  is  poor  and  the  discrepancy  increases  as  the  normal  size  becomes  smaller. 
This  poor  grid  convergence  probably  is  because  this  spanwise  mesh  size  is  too  large  for  the  elongated  turbulence 
structures  very  near  the  wall  in  this  subsonic  boundary- layer  flow.  Grid  sizes  should  be  designed  according  to  the 
turbulent  structures. 


IV.  Boundary-Layer  Flows  in  a  C-D  Nozzle 

JENRE  has  been  used  frequently  to  study  supersonic  jet  noise  at  various  jet  operating  conditions  and  nozzle 
geometries.  For  example,  the  noise  generation  in  imperfectly  expanded  jets  [8],  the  noise  propagation  in  the 
presence  of  a  deck  and  a  jet  blast  deflector  [9],  and  the  assessment  of  noise  reduction  using  chevrons  [10].  The  wall- 
model  approach  is  thus  tested  in  boundary  layers  inside  the  C-D  nozzle  used  in  those  studies.  It  should  be  mentioned 
that  there  are  non-zero  pressure  gradients  inside  the  nozzle,  so  the  flow  is  not  truly  equilibrium.  But  it  may  be 
argued  that  the  non-equilibrium  effect  is  accounted  for  by  the  velocity  of  the  first  node  above  the  wall.  The 
equilibrium  wall-model  method  has  been  used  previously  by  some  researchers  to  simulate  the  boundary- layer  effect 
in  jet  noise  simulations,  such  as  the  work  by  Aikens,  et  al.  [15]  and  Deniau,  et  al.  [16].  Both  studies  use  finite- 
difference  solvers  and  both  observed  the  numerical  issue  associated  with  using  the  first  point  above  the  wall.  On  the 
other  hand,  Guillaume  et  al  used  a  hybrid  method  to  implement  the  boundary- layer  effect  in  their  jet  noise 
simulations  using  a  finite-volume  solver  [17].  As  we  have  mentioned  in  the  Introduction  section,  the  wall-model 
implementation  in  a  finite-element  solver  for  large-eddy  simulations  has  rarely  been  reported.  To  our  best 
knowledge,  this  is  the  first  time  the  application  of  a  wall-model  method  to  a  finite-element  solver  being  presented  to 
the  jet  noise  simulation  community. 
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The  nozzle  geometry  is  shown  in  Figure  12,  which  is  representative  of  practical  supersonic  military  engine  nozzles. 
The  design  Mach  number  Md  of  this  C-D  nozzle  is  1.5,  and  the  fully  expanded  pressure  ratio  is  3.7.  The  nozzle 
pressure  ratio  ( NPR )  is  4.0,  which  produces  an  underexpanded  nozzle-exit  condition.  The  Reynolds  number  based 
on  the  nozzle  exit  condition  is  slightly  above  four  million.  The  boundary- layer  thickness  (5)  near  the  nozzle  lip  is 
roughly  1.2%  of  the  nozzle-exit  diameter  (D).  Because  the  Reynolds  number  is  high  and  the  boundary  layer  is  very 
thin,  it  is  not  practical  to  place  many  cells  inside  the  boundary  layer.  Thus,  tetrahedral  meshes  with  large  aspect 
ratios  in  both  streamwise  and  spanwise  directions  are  used  to  reduce  the  computational  cost.  The  cell  size  in  the 
principal  flow  direction  is  0.35  and  0.335  in  the  circumstantial  direction.  They  are  indeed  very  large  cell  sizes,  much 
larger  than  those  reported  in  other  wall-model  implementations  [6] [15] [16].  Three  cell  sizes  in  the  normal  direction 
are  tested  and  they  are  shown  in  Table  2.  The  growth  rate  of  the  cell  size  is  set  to  1.05  in  the  two  meshes  where  the 
normal  grid  size  near  the  wall  is  smaller  than  the  streamwise  cell  size.  Because  even  the  finest  grid  resolution  only 
has  11-12  cells  inside  the  boundary  layer,  all  the  three  normal  grid  resolutions  are  coarse  in  terms  of  the  boundary- 
layer  simulation.  Surface  roughness  is  implemented  in  the  subsonic  region  upstream  of  the  C-D  section  to  trip  the 
boundary  layer.  The  roughness  height  is  1%  of  the  nozzle  exit  diameter,  which  is  comparable  to  the  boundary- layer 
thickness  near  the  nozzle  lip.  Figure  13  shows  the  instantaneous  vorticity  contours  inside  and  around  the  nozzle. 
Vortex  structures  are  generated  by  the  surface  roughness,  and  they  are  stretched  when  they  flow  through  the 
diverging  section  of  the  C-D  nozzle  due  to  the  favorable  pressure  gradient.  Figure  14  shows  the  velocity  profiles  at 
two  axial  locations  upstream  of  the  nozzle  lip.  There  is  a  very  good  grid  convergence  when  the  wall  model  is 
applied.  The  velocity  profiles  align  well  with  the  logarithmic  distribution  given  in  Eq.  (2).  However,  the  unresolved 
boundary  layers  where  the  wall  model  is  not  implemented  present  much  larger  velocity  deficits,  which,  as  expected, 
decreases  as  the  normal  grid  size  decreases.  It  appears  that  the  mesh  type  that  has  both  large  streamwise  and 
spanwise  aspect  ratios  performs  well  in  this  supersonic  nozzle  flow.  This  may  indicate  that  turbulence  structures  in 
the  outer  layer  given  in  these  simulations  have  spanwise  and  streamwise  lengths  larger  than  the  specified  spanwise 
and  streamwise  cell  sizes.  Figure  15  shows  turbulence  intensity  profiles  of  the  axial  velocity  slightly  downstream  of 
the  nozzle  exit.  Those  results  are  interpolated  from  the  data  of  unstructured  grid  nodes.  It  can  be  seen  that  the 
profiles  predicted  by  the  three  grid  resolutions  are  similar,  except  the  peak  location  in  the  radial  direction.  The  finer 
grid  resolution  predicts  a  peak  radial  location  further  from  the  jet  center.  This  is  consistent  with  the  observation  that 
the  peak  turbulence  intensity  shifts  towards  the  wall  as  the  grid  resolution  increases.  It  should  also  be  mentioned 
that  the  current  wall-model  implementation  observes  a  grid-to-grid  oscillation  in  the  axial  direction  in  supersonic 
flows  when  the  structured-like  mesh  is  used.  The  oscillation  magnitude  increases  as  the  grid  aspect  ratio  increases, 
and  a  spatial  averaging  is  needed  to  eliminate  this  numerical  oscillation.  But  this  numerical  phenomenon  is  not 
present  when  an  isotropic  mesh  is  used.  This  probably  is  because  the  first  node  adjacent  to  the  wall  is  not  directly 
above  the  geometric  center  of  the  cell  face  on  the  wall  in  a  structured-like  mesh,  but  this  is  more  likely  so  in  an 
isotropic  mesh. 


V.  Conclusions 

The  equilibrium  wall  model  is  implemented  in  our  in-house  finite  element  flow  solver  JENRE  to  simulate  the 
boundary-layer  effect  in  flows  at  high  Reynolds  numbers.  The  Crocco-Busemann  relation  is  used  to  account  for  the 
compressibility.  The  no-slip  adiabatic  boundary  condition  is  applied  to  the  inviscid  and  viscous  fluxes  on  the  wall  to 
satisfy  the  surface  physical  condition.  However,  since  the  linear  basis  function  is  used  in  JENRE,  a  zero-surface- 
velocity  would  greatly  underestimate  the  volume  and  surface  integrals  in  the  first  cell  adjacent  to  the  wall  comparing 
to  those  derived  from  the  logarithmic  velocity  profile  and  the  Crocco-Busemann  relation.  Thus,  the  surface 
tangential  velocity  is  not  forced  to  zero  in  those  integrations.  This  surface  tangential  velocity  is  updated  along  with 
the  flow  field  at  interior  points.  This  wall-model  implementation  is  validated  in  a  subsonic  boundary-layer  flow  at 
Moo  =  0.9  over  a  flat  plate  and  a  supersonic  flow  in  a  C-D  nozzle,  which  is  used  frequently  in  our  previous  jet  noise 
simulations.  The  inflow  turbulences  are  generated  by  surface  roughness  in  both  cases.  It  is  found  that  this 
implementation  presents  an  excellent  grid  convergence.  It  does  not  present  the  numerical  issue  associated  with  using 
the  first-point  as  reported  in  other  wall-model  implementations.  Grid  sizes  much  larger  than  those  recommended  in 
other  wall-model  implementations  can  be  used.  Skin  frictions,  velocity  profiles  and  turbulence  quantities  agree  well 
with  available  experimental  data  and  theoretical  models.  In  addition,  there  is  a  very  good  agreement  between  the 
shape-factor  predictions  and  those  derived  from  the  equilibrium  turbulent  boundary-layer  velocity  profile. 
Furthermore,  there  is  a  good  agreement  between  the  predicted  volume  and  surface  integrals  in  a  cell  adjacent  to  the 
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wall  based  on  a  non-zero  tangential  surface  velocity  and  those  derived  from  the  equilibrium  turbulent  boundary- 
layer  velocity  profile  and  the  density  profile  based  on  the  Crocco-Busemann  relation. 
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Figure  1.  Velocity  profile  in  an  equilibrium  turbulent  boundary  layer.  The  definition  of  u+,y+,  k  and  B  can  be  found 
in  the  section  II. 


7 


Table  1.  Mesh  parameters  at  Ree  =  31000  (using  cubic  cells) 
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Figure  3.  An  instantaneous  temperature  distribution  (using  Mesh_II  in  Table  1)  of  a  boundary-layer  flow  at  = 
0.9  on  a  flat  plate. 


Figure  4.  Axial  evolvement  of  the  skin-friction  in  a  boundary-layer  flow  at  =  0.9  on  a  flat  plate.  The  description 
of  simulation  data  can  be  found  in  Table  1 .  The  skin  friction  and  the  Reynolds  number  ReQ  are  equivalent 

incompressible  values  estimated  by  the  van  Driest  II  transformation  [13]  shown  in  Eq.  (4).  Dashed  line:  K'arm'an- 
Schoenherr  incompressible  empirical  correlation  [18].  Symbols:  incompressible  experimental  data  from  De  Graaff 
and  Eaton  [19]. 
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Figure  5.  Streamwise  velocity  profiles.  Predictions  using  the  three  coarser  meshes  are  at  the  equivalent 
incompressible  Reeinc  =  31000.  The  result  of  the  finest  mesh  is  at  Reeinc  =  20000.  The  description  of  simulation 
data  can  be  found  in  Table  1.  Dashed  line:  logarithmic  velocity  profile  shown  in  Eq.  (2).  Symbols:  incompressible 
experimental  data  from  De  Graaff  and  Eaton  at  Ree  =  31000  [19]. 


Figure  6.  Comparison  of  turbulence  intensities  and  the  Reynolds  shear  stress  at  Red>inc  =  31000.  All  quantities  are  in 
wall  units.  The  description  of  simulation  data  can  be  found  in  Table  1.  Symbols:  incompressible  experimental  data 
from  De  Graaff  and  Eaton  at  Ree  =  31000  [19]. 
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Figure  7.  Axial  distributions  of  the  shape  factor  in  a  boundary-layer  flow  at  Mm  =  0.9  on  a  flat  plate.  The 
description  of  simulation  data  can  be  found  in  Table  1.  (a)  In  the  first  cell  above  the  wall,  the  velocity  and  density 
are  assumed  to  vary  linearly  with  the  vertical  distance  from  the  wall,  (b)  In  the  first  cell  above  the  wall,  the  velocity 
uses  the  equilibrium  boundary-layer  velocity  profile  shown  in  Figure  1  and  the  density  profile  follows  the  Crocco- 
Busemann  relation. 
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Figure  8.  Ratios  between  the  predicted  momentum  integrals  (J  pu  dy )  in  the  first  cell  above  the  wall  and  those 
derived  from  the  equilibrium  boundary- layer  velocity  profile  (Figure  1)  and  the  density  profile  based  on  the  Crocco- 
Busemann  relation.  The  solid  lines:  using  a  non-zero  surface  velocity.  The  dashed  lines:  using  a  zero-surface 
condition. 


Figure  9.  Ratios  between  the  predicted  kinetic  energy  integrals  (J  pu2  dy)  in  the  first  cell  above  the  wall  and  those 
derived  from  the  equilibrium  boundary- layer  velocity  profile  (Figure  1)  and  the  density  profile  based  on  the  Crocco- 
Busemann  relation.  The  solid  lines:  using  a  non-zero  surface  velocity.  The  dashed  lines:  using  a  no-slip  condition. 
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Figure  10.  Comparison  of  streamwise  velocity  profiles  at  Reeinc  =  31000.  The  cell  size  in  the  streamwise  direction 
is  fixed  at  Ax+  =  370.  The  normal  and  the  spanwise  grid  sizes  (Ay  =  Az]  vary,  and  there  are  three  aspect  ratios 
(Ax/Ay).  The  description  of  simulation  data  can  be  found  in  Table  1.  Dashed  line:  logarithmic  velocity  profile 
shown  in  Eq.  (2).  Symbols:  incompressible  experimental  data  from  De  Graaff  and  Eaton  at  Ree  =  31000  [19]. 


Figure  11.  Comparison  of  streamwise  velocity  profiles  at  Ree>inc  =  31000.  The  cell  sizes  in  both  the  streamwise  and 
spanwise  directions  are  fixed  at  Ax+  =Az+  =  370.  The  normal  cell  size  varies,  and  there  are  three  aspect  ratios 
(Ax/Ay).  The  description  of  simulation  data  can  be  found  in  Table  1.  Dashed  line:  logarithmic  velocity  profile 
shown  in  Eq.  (2).  Symbols:  incompressible  experimental  data  from  De  Graaff  and  Eaton  at  Ree  =  31000  [19]. 
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Figure  12.  The  converging  and  diverging  nozzle  geometry,  the  wall  roughness  shape  and  the  boundary-layer  mesh 
near  the  nozzle  surfaces. 


Table  2.  Mesh  parameters  of  the  boundary-layer  mesh  slightly  upstream  of  the  nozzle  lip. 
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Figure  13.  An  instantaneous  vorticity  distribution  cd/(Uj-D). 
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Figure  14.  Axial  velocity  profiles  at  two  axial  locations  upstream  of  the  nozzle  exit,  “h”  is  the  vertical  cell  size  Ay  of 
the  first  cell  adjacent  to  the  wall.  5  is  the  boundary-layer  thickness,  (a)-(b)  linear  profiles.  Solid  lines  are  predictions 
using  the  wall-model  method.  Dashed  lines  are  velocities  of  unresolved  boundary  layers,  where  the  wall  model  is 
not  implemented,  (c)-(d)  logarithmic  profiles.  The  axial  velocities  are  the  van  Deiest  effective  velocities  [13].  The 
dash-dotted  line  is  the  logarithmic  velocity  profile  defined  in  Eq.  (2). 


Figure  15.  Radial  turbulence  intensity  profiles  slightly  downstream  of  the  nozzle  exit  at  x  =  0.1D.  “h”  is  the  vertical 
cell  size  Ay  of  the  first  cell  adjacent  to  the  wall.  5  is  the  boundary- layer  thickness. 
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